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ABSTRACT 

A numerical model of HgCdTe solidification was implemented using finite the element code FIDAP. Model 
verification was done using both experimental data and numerical test problems. The model was used to 
evaluate possible effects of double-diffusion convection in molten material, and microgravity level on 
concentration distribution in the solidified HgCdTe. Particular attention was paid to incorporation of 
HgCdTe phase diagram. It was found, that below a critical microgravity amplitude, the maximum convective 
velocity in the melt appears virtually independent on the microgravity vector orientation. Good agreement 
between predicted interface shape and an interface obtained experimentally by quenching was achieved. The 
results of numerical modeling are presented in the form of video film. 

Keywords: Binary alloy, HgCdTe, bulk crystal growth, directional solidification, microgravity, numerical 
modelling, FIDAP. 

1. INTRODUCTION 

Solidification of HgCdTe was the subject of recent numerical modeling 14 . In these papers particular 
attention was paid to proper modeling of melt convection and solidification to predict the crystal/melt 
interface shape and crystal radial and axial inhomogenity. In order to achieve this, HgCdTe was considered 
as binary alloy with melting temperature available from phase diagram. The present study is a continuation 
of these efforts with the emphasis on investigation of melt convection on microgravity level, double- 
diffusion and interface shape. It is useful to clarify from the beginning, that information about actual crystal 
/ melt interface during growth process is available a posteriori, and can be obtained from radial distribution 
of mole concentration of CdTe and from visible change in crystal morphology, that occurs if the growth 
process is interrupted by fast sample removal from the furnace (quench). No direct information about 
convection velocity is available. In the present work double diffusion convection in HgCdTe melt and 
crystal-melt interface curvature prediction during directional solidification are addressed numerically using 
the finite element approach. * 

It appears that both global modeling of heat transfer in the furnace along with a detailed solidification model 
are required to fulfill proper modeling of crystal growth as well as furnace. A global model of AADSF 
furnace , that includes radiative, conductive and convective heat transfer was used to obtain the appropriate 
temperature boundary conditions for detailed solidification model. Alternatively experimental temperature 
measurements were used to generate temperature boundary conditions for the numerical model. The 
numerical model of convection and solidification of binary alloy is based on the general equations of heat 
and mass transfer in two-dimensional region. Minimal code development was undertaken, as commercial 
finite element code FIDAP 6 was used. Two complimentary solutions where obtained: one for melt region 
with a fixed boundary and other for melt region with an unknown boundary. In the former case, numerically 
converged solutions where typically easily obtained with standard FIDAP settings, while in the latter case a 
special choice of relaxation factors for each degree of freedom was required. The simplified model with 
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fixed interface shape was especially useful for parametric study of convection, as the steady state solution 
was available for a broad range of parameters. 

In the complete numerical model with unknown boundary the following thermophysical properties of 
HgCdTe was taken into account: the pseudobinary HgTe-CdTe phase diagram 7 , the compositional 
dependence of the interface segregation coefficient, an assumed diffusion coefficient for mass transport in 
HgCdTe, the temperature and composition dependence of the melt density (to calculate buoyancy force) and 
the melt viscosity and its temperature dependence 8 . The composition dependence of diffusion coefficient 
was not considered, although it might be important for double-diffusion convection. In the simplified model 
with the fixed boundary melting temperature was fixed instead of being calculated from phase diagram. 

The basic physics that controls solidification front shape formation is as follows. Due to the large separation 
between the liquidus and the solidus, there is a wide variation of composition through the boundary layer 
adjacent to the interface. Cooling from below and the segregation at the solid-melt interface of the heavier 
HgTe-rich solute each tend to suppress convection. However, the thermal conductivity of the solid is a factor 
of 8 times smaller then the melt (see table 1) and there is a resulting thermal short circuit through the 
ampoule (fused silica); this results in curved isotherms in this region as the heat cannot travel parallel to the 
axis within the center of the sample. The interface is therefore curved and, due to the thermal and solutal 
convection is neither an isotherm nor a line of isoconcentration. Double diffusion also contributes to the 
mass transport. In the present study, the numerical model is used to consider both these source of mass 
transport. 

2. MATHEMATICAL MODEL AND NUMERICAL PROCEDURE 

The main features of the model are double-diffusion convection in the melt and solidification. The set of 
steady state Navier-Stokes equations in Boussinesq approximation along with energy and solutal balance 
equations are used 6 . In the following u is melt convection velocity, R is sample pulling velocity (assumed at 
the steady state to be equal to crystal growth velocity), p - pressure, p - density, g -microgravity acceleration, 
P T and Pc - volumetric expansion coefficients due temperature and concentration change respectively, T 0 
and Co -reference temperature and concentration, p - dynamic viscosity, D - mass diffusion coefficient, k - 
thermal conductivity, Cp - specific heat at constant pressure. Furthermore indexes “solid" and “melt" are 
used to refer to solidified and molten parts of the sample. 

Momentum conservation equation (1) 

p(mVm) = -V/? + jj.A« + pg[i -p 7 (r-r 0 ) + p c (c-c 0 )] 


Mass conservation equation (Boussinesq approximation) 


( 2 ) 


Vw = 0 


Energy conservation equation 


(3) 


pC p {uVT) = W(kVT) 

Species conservation equation (4) 


wVC = DAC 


242 



The boundary conditions incorporate major physical assumptions in order to describe solidification of an 
unknown interface between molten and solidified HgCdTe. These phase change boundary conditions are as 
follows. 

The melting temperature is determined by (5) 

T melt ~ T solid ~ T 0^ C 0^ + m ^ C melt 
where (6) 

dT_ 

dC 

The heat balance on the interface, including latent heat release (L) is (7) 

k melt^ T melt ' n ~ k solid^ T solid ' n = P solid L [^ ~ jj) ' n 

Here t stands for pseudo-time, i.e. artificial melting time, used in the model to allow interface adjustment in 
accordance with temperature field and phase diagram. Notation n is unit vector normal to the interface. 

The balance of the mass flux across interface is (8) 

= p 

The no-slip boundary condition at the interface is (9) 

(m - R) x n — 0 

The mass balance for solute transport across the interface is (10) 



^melt 


P solid melt 

= c 


D V C melt ' " “ D solid^ C solid ' " 


„ dS\ » r 

solid{ K ~dt)' n melt 


dS\ * 
U ~di)' n 


The chemical equilibrium on the interface is 


C — k (C)C 

solid segregation v f melt 


(ID 

( 12 ) 


where 

(r\ - Solidus(C 7) 
K segregation t > ~ Liquidus(C) 


(13) 


To perform actual calculations initial temperature and solutal field, convection velocities and interface 
position are required. In order to obtain a steady state solution of the problem (1*4) with boundary 
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conditions (5-11) a three step numerical procedure was used. On the first computational step the position of 
interface was flat and fixed, and so variable S = 0. Initial zero velocity field and constant initial distribution 
of temperature and concentration were used to obtain the solution of double-diffusion convection problem in 
two-dimensional area with fixed boundaries, that satisfies temperature boundary conditions obtained from 
global numerical model of AADSF furnace or using available measured temperature profile. This solution 
was used as initial conditions during the second step, when variable S was released and the interface shape 
was obtained with proper adjustment of the computational mesh. A final third step was required to obtain a 
convergent solution for all variables involved. Two-dimensional 9-node quadrilateral finite elements were 
implemented. Due to strong nonlinearity of the problem, a relaxation factor value as high as 0.9 was used for 
surface elevation variable S. A first-order upwinding numerical scheme was used. This ensures suppression 
of “wiggles”, a typical problem caused by inadequate mesh density to resolve strong gradients of 
concentration. Further details on finite element implementation of the model can be found elsewhere . For 
parametric studies a simplified model was used. The major simplification was the elimination of 
independent variable S with appropriate boundary conditions on the fixed solid/ melt interface. 

3, NUMERICAL TEST PROBLEM: DOUBLE-DIFFUSION CONVECTION 

In order to verify the numerical approach and choice of numerical parameters, such as mesh size and 
relaxation factors, the following problem was considered 9 . Double-diffusion convection governed by set of 
equations 1-4 occurs in a two-dimensional enclosure which is heated and salt rich at the bottom. This is a a 
case of opposing or counteracting flow, so that as the heat impact tends to start convection, the salt 
introduction from below tends to suppress it. It is also known as a “diffusive” regime, as the component with 
a higher diffusivity (heat) has a destabilizing effect, while the component with a lower diffusivity (salt) has a 
stabilizing effect on the local density gradient. Instead of the complicated boundary conditions 5-11, 
representing crystal growth in the region with the unknown boundary, a straightforward set of boundary 
conditions apply. Constant temperature and salinity on the top and the bottom of enclosure 
are determined by the 0^) 
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while the isolated vertical boundaries of enclosure are 

x = 1, 0 

dy 

(15) 


x = 1, 0 

~ = ° 
dy 



In this case the following non dimensional parameters, which control fluid behavior can be 
introduced. Thermal and solutal Rayleigh numbers, based on enclosure height H and top to bottom 
differences of temperature and salt are as follows 


(rj-Ty-Pr *-# 3 p (c,-c 0 )P c s-h 3 
Ra r = Ra c = 
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ermal Rayleigh numberResults of test problem appears 







Here k x and k c - are thermal and mass diffusivity, respectively and v is kinematic viscosity. The complete 
set of nondimesional parameters for this problem include also Prandtl and Lewis numbers. Geometry is 
characterized by aspect ratio. 



To monitor convection influence on heat transfer we will use an overall Nusselt number at vertical walls 

r h .. , XT H (dT\ 

Nu = L Nu^dy where Nu = — 3 - 

}0 y J y T l -T 0 \dyjwall 

Similarly the impact of convection on solutal flux is characterized by overall Sherwood number 

Sh = f”sh y dy .Here Sh y = 

In order to perform comparison with published 9 results the nondimensional parameters where fixed, as 
follows: Ra c =2000, Pr=l, Le=10, H/L=l. Thermal Rayleigh number was varied in wide range from 2000 to 
18000 in order to determine the critical thermal Rayleigh number, when convection occurs. Results 
presented on fig 1, indicate good agreement with data from 9 . 

4. RESULTS AND DISCUSSION 
Two complimentary solutions are presented here. 

The first one is based on the complete model described in section 2. In this case typical distribution of 
temperature field along with pattern of melt flow is presented on fig2a. For comparison experimental 
photograph featuring “quenched” interface is presented on fig 2b. Good agreement in shape is mainly due 
proper choice of thermal conductivities of solidified and molten material as well as silica ampoule thermal 
conductivity and wall thickness. The relevant material and system parameters are presented in the tables 1-3. 
This is the case, the strict vertical orientation of the sample manifests itself in a symmetrical pattern of the 
flow. In other simulated cases deviation from the vertical of only 1 degree cases loss of symmetry and 
formation of one cell pattern of convection. These cases are featured in the video film. The shape of interface 
closely fits and is almost the same for various levels of gravity. This strongly suggests that convection in the 
melt virtually has no effect on solidification interface curvature in this particular case. Further efforts are 
required to investigate this finding in more detail. It should be noted, that simulation of this problem using 
non-dimensional parameters, like Rayleigh number, can be confusing in this case, as the driving force for 
convection is not a vertical or a simple radial gradient, but rather a radial temperature gradient created by the 
curvature of the solidification front. Otherwise the cooling of the melt from below and rejection of the 
heavier component (HgTe) would suppress convection. In this case both the shape of and concentration 
distribution along the interface are of major interest. Radial distribution of concentration and temperature are 
presented on fig 3. The radial variation in concentration is only in fair agreement with experimental 
observation. It appears, that the calculated temperature radial variation is about 7°C less than expected from 
the experimental observed variation of 20 °C. 

The second set of solutions is based on simplified model with fixed solid / melt interface. Parametric study 
of maximum convection velocity depending on microgravity conditions is presented on fig 4a and 4b. 
Geometry of the sample and temperature boundary conditions was taken from experiment SLM13. While 
material properties was the same as in former case, the system parameters was different, as shown in table 3. 
Figure 4a features convection velocity dependence on microgravity amplitude. Figure 4b presents the same 
data depending on gravity vector angle to sample axis. It is clear, that the maximum velocity of convection 
occurs, when ampoule is inclined at 90°, while if ampoule is strictly vertical convection is less by factor of 
50%. These data strongly suggest, that in order to reduce convection in the melt, vertical orientation of the 






Table 1. HgCdTe thermophysical property data 


Thermal conductivity 
(solid) 

0.291 

W/°Km 



Thermal conductivity 
(melt) 

1.96 

W/°Km 



Density (solid) 

7630 

kg / m 3 



Density (melt) 

7550 

kg / m 3 



Specific heat (solid) 

177 

J / °K kg 



Specific heat (melt) 

257 

J / °K kg 



Kinematic viscosity 

1.08 x 10' 7 

m 2 / s 



Latent heat 

130000 

J/kg 



Thermal expansion 
coefficient 

1.5 x 10' 4 

° K -i 



Solutal expansion 
coefficient 

0.3 

°K-i 



Diffusion coefficient of 
CdTe 

5.5 x 10' 9 

m 2 /s 



Melting temperature 
(liquidus) 

944.09 + 671.46 x C - 375.07 x C 2 + 
110.54 xC 3 

°K 



Segregation coefficient 

4.7 - 15.9 x C + 39.3 x C 2 -59.1 x C 3 + 





46.7 xC 4 - 14.7 xC 5 




Table 2. Ampoule (Quartz) thermophysical property data 




Thermal conductivity 

2.5 

W/°Km 


Density 

2200 

kg / m 3 


Specific heat 

1050 

J / °K kg 


Table 3. System data 





Case 

SL-413 

MCT-16Q 


i 

Sample radius 

3.9 

3.9 

mm 

Sample length 

110 

138 

mm 

Ampoule outer radius 

5.9 

5.9 

mm 

Pulling velocity 

2.0x 10 ' 7 

2.0xl0' 7 

m/s 

Cold heater temperature 

312 

360 

°c 

Hot heater temperature 

827 

846 

°c 

CdTe initial concentration 

0.152 

0.2 


j 





Concentration 


0.15 


Figure 3. 


Raduis R/Ro 


Concentration and temperature distribution on the solidification interface. Crosses 
distribution (axis on the left). Circles - temperature distribution (axis on the right). 


concentration 


Temperature C 






Figure 4. 

a) Maximum convection velocity as function of microgravity acceleration amplitude. 

b) Maximum convection velocity as function of gravity vector angle to sample axis. Angle of 90 degree 
corresponds to the strict vertical orientation of the sample. 
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ampoule is of utmost importance. This is valid till gravity level of order 10 6 of gravity on Earth. However at 
very low gravity, of order 10~ 7 and 10' 8 of earth gravity the tilt of the ampoule has much less effect. 

5. CONCLUSIONS 

Prediction of solidification interface shape and flow pattern in molten HgCdTe for various gravity conditions 
are available using a 2-D model of melt convection and solidification. This model utilizes temperature 
boundary conditions from results of modeling of heat transfer in AADSF furnace, and incorporates all basic 
modes of heat transfer, such as conduction, convection and radiation. Alternatively, a measured experimental 
temperature profile can be used. 

It appears that crystal/melt interface shape, as well as concentration distribution on the interface depend 
primarily on heat flow and concentration distribution in close vicinity of interface. The models ability to 
predict correctly the interface shape and concentration distribution heavily depends on known material 
properties, such as conductivity and segregation coefficient. 

Maximum convection velocity in the ampoule is extremely sensitive to gravity vector orientation and can be 
reduced at least by factor of 50% if vertical orientation of ampoule is strictly observed. This conclusion is 
valid to microgravity level of about 10 7 compared to earth gravity. If microgravity is lower than this level 
convection in the melt is not anymore sensitive to the ampoule orientation. 
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